library(tidyverse)
library(data.table)
library(ggh4x)
library(lme4)
library(car)
library(ggeffects)
library(doParallel)

cores <- getOption("mc.cores", detectCores()) # for parallel computation
cl <- makeCluster(cores)
registerDoParallel(cl)

data loading

files <- dir("exp1_data_eyetracking_individual_saccades", ".*csv$", full.names = TRUE) 
dat <- c()

for (i in 1:length(files)) {
    d <- fread(files[i], header = TRUE)
    d$id <- i  # numeric variable for subject id
    dat <- rbind(dat, d)
}

dat <- subset(dat, dat$event == "fixation" & dat$Conf != 0)
dat <- mutate(dat, target = apply(dat, 1, function(x){names(d)[13 + which.max(x[14:16])]}))
dat <- mutate(dat, dud = apply(dat, 1, function(x){names(d)[13 + which.min(x[14:16])]}))
dat$Condition <- as.factor(dat$Condition)
dat$target <- fct_recode(dat$target, up = "UVal", left = "LVal", right = "RVal")
dat$dud <- fct_recode(dat$dud, up = "UVal", left = "LVal", right = "RVal")
dat <- mutate(dat, fix = ifelse(item == target, "target", 
                         ifelse(item == dud, "dud",
                         ifelse(item == "noFix", "noFix", "distractor"))))
dat <- mutate(dat, choice = ifelse(ChosenITM == target, "correct", 
                            ifelse(ChosenITM == dud, "dud", "distractor")))

subject-wise fixation plot

plot1 <- foreach(i = 1:length(files), .packages = c("tidyverse", "ggh4x")) %dopar% {
    subset(dat, dat$id == i) %>%
    ggplot() + geom_point(aes(x = x, y = y, size = dur, color = Condition), alpha = 0.3) + 
        facet_nested(. ~ target + dud) + ggtitle(i) -> p
    print(p)
}
plot1
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

Position-based fixation frequency

dat %>%
    group_by(Condition, target, dud, item, id) %>%
    summarise(n = n()) -> freq
## `summarise()` has grouped output by 'Condition', 'target', 'dud', 'item'. You
## can override using the `.groups` argument.
plot2 <- foreach(i = unique(freq$Condition), .packages = c("tidyverse", "ggh4x")) %dopar% {
    p <- ggplot(subset(freq, freq$Condition == i)) + geom_violin(aes(x = item, y = n, color = item)) + geom_point(aes(x = item, y = n)) + 
        facet_nested(. ~ target + dud) + 
        scale_x_discrete(limits = c("up", "left", "right", "noFix")) +
        scale_color_discrete(limits = c("up", "left", "right", "noFix")) + ggtitle(i)
    print(p)
}
plot2
## [[1]]
## Warning: Groups with fewer than two data points have been dropped.

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

Stimulus-based fixation frequency

dat %>%
    group_by(Condition, fix, id) %>%
    summarise(n = n()) -> fixation
## `summarise()` has grouped output by 'Condition', 'fix'. You can override using
## the `.groups` argument.
fixation$fix <- as.factor(fixation$fix)
fixation$fix <- factor(fixation$fix, levels = c("target", "distractor", "dud", "noFix"))

ggplot(fixation, aes(x = fix, y = n, color = Condition)) + geom_violin() + 
    geom_point(position = position_dodge(width = .9)) +
    stat_summary(fun.y = "mean", geom = "crossbar", position = position_dodge(width = .9)) 
## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.

f1 <- lmer(n ~ fix * Condition + (1 + fix + Condition | id), data = fixation)
## boundary (singular) fit: see help('isSingular')
summary(f1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: n ~ fix * Condition + (1 + fix + Condition | id)
##    Data: fixation
## 
## REML criterion at convergence: 2034.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7499 -0.4896 -0.0257  0.4725  4.2569 
## 
## Random effects:
##  Groups   Name          Variance  Std.Dev. Corr                               
##  id       (Intercept)    7389.117  85.960                                     
##           fixdistractor     1.926   1.388   1.00                              
##           fixdud         3687.390  60.724  -0.97 -0.97                        
##           fixnoFix      13342.569 115.510  -0.92 -0.92  0.86                  
##           Condition0.3     40.107   6.333   0.49  0.49 -0.42 -0.76            
##           Condition0.5     92.452   9.615   0.91  0.91 -0.79 -0.94  0.62      
##           Condition0.7     64.211   8.013   0.88  0.88 -0.83 -0.96  0.63  0.90
##           Condition0.85   226.407  15.047   0.75  0.75 -0.72 -0.91  0.75  0.76
##           Condition0.95   282.650  16.812   0.76  0.76 -0.73 -0.88  0.60  0.77
##  Residual                 369.516  19.223                                     
##             
##             
##             
##             
##             
##             
##             
##             
##   0.95      
##   0.96  0.98
##             
## Number of obs: 239, groups:  id, 10
## 
## Fixed effects:
##                             Estimate Std. Error t value
## (Intercept)                  393.700     27.854  14.134
## fixdistractor                -22.400      8.608  -2.602
## fixdud                      -379.911     21.150 -17.963
## fixnoFix                      -5.900     37.525  -0.157
## Condition0.3                 -16.700      8.827  -1.892
## Condition0.5                 -26.500      9.119  -2.906
## Condition0.7                 -41.600      8.962  -4.642
## Condition0.85                -34.600      9.826  -3.521
## Condition0.95                -33.700     10.108  -3.334
## fixdistractor:Condition0.3     6.600     12.158   0.543
## fixdud:Condition0.3           60.811     12.349   4.925
## fixnoFix:Condition0.3         27.200     12.158   2.237
## fixdistractor:Condition0.5     2.800     12.158   0.230
## fixdud:Condition0.5          110.111     12.349   8.917
## fixnoFix:Condition0.5         39.800     12.158   3.274
## fixdistractor:Condition0.7    -0.100     12.158  -0.008
## fixdud:Condition0.7          200.811     12.349  16.262
## fixnoFix:Condition0.7         54.700     12.158   4.499
## fixdistractor:Condition0.85  -15.200     12.158  -1.250
## fixdud:Condition0.85         254.411     12.349  20.602
## fixnoFix:Condition0.85        46.000     12.158   3.784
## fixdistractor:Condition0.95  -30.800     12.158  -2.533
## fixdud:Condition0.95         302.111     12.349  24.465
## fixnoFix:Condition0.95        40.100     12.158   3.298
## 
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(f1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: n
##                  Chisq Df Pr(>Chisq)    
## fix            569.807  3  < 2.2e-16 ***
## Condition       50.098  5  1.323e-09 ***
## fix:Condition 1390.787 15  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(f1, terms = c("fix", "Condition"), type = "fe"))

Stimulus-based fixation frequency (choice considered)

dat %>%
    group_by(Condition, fix, choice, id) %>%
    summarise(n = n()) -> fixation2
## `summarise()` has grouped output by 'Condition', 'fix', 'choice'. You can
## override using the `.groups` argument.
fixation2$fix <- as.factor(fixation2$fix)
fixation2$fix <- factor(fixation2$fix, levels = c("target", "distractor", "dud", "noFix"))

# correct choice
subset(fixation2, fixation2$choice == "correct") %>%
    ungroup() %>% complete(Condition, fix, choice) %>%
    ggplot() + geom_point(position = position_dodge(width = .8), aes(x = fix, y = n, color = Condition)) +
    stat_summary(fun.y = "mean", geom = "crossbar", position = position_dodge(width = .8), 
                 mapping = aes(x = fix, y = n, color = Condition)) + ggtitle("fixation frequency in correct trials")

# distractor choice
subset(fixation2, fixation2$choice == "distractor") %>%
    ungroup() %>% complete(Condition, fix, choice) %>%
    ggplot() + geom_point(position = position_dodge(width = .8), aes(x = fix, y = n, color = Condition)) +
    stat_summary(fun.y = "mean", geom = "crossbar", position = position_dodge(width = .8), 
                 mapping = aes(x = fix, y = n, color = Condition)) + ggtitle("fixation frequency in distractor-chosen trials")

# dud choice
subset(fixation2, fixation2$choice == "dud") %>%
    ungroup() %>% complete(Condition, fix, choice) %>%
    ggplot() + geom_point(position = position_dodge(width = .8), aes(x = fix, y = n, color = Condition)) +
    stat_summary(fun.y = "mean", geom = "crossbar", position = position_dodge(width = .8), 
                 mapping = aes(x = fix, y = n, color = Condition)) + ggtitle("fixation frequency in dud-chosen trials")
## Warning: Removed 6 rows containing non-finite values (`stat_summary()`).
## Warning: Removed 6 rows containing missing values (`geom_point()`).

Stimulus-based fixation frequency (confidence considered)

dat %>%
    group_by(Condition, fix, choice, Conf, id) %>% 
    summarise(n = n()) %>% ungroup() %>% complete(Condition, fix, Conf) -> fixation3
## `summarise()` has grouped output by 'Condition', 'fix', 'choice', 'Conf'. You
## can override using the `.groups` argument.
fixation3$fix <- as.factor(fixation3$fix)
fixation3$fix <- factor(fixation3$fix, levels = c("target", "distractor", "dud", "noFix"))

# correct choice
subset(fixation3, fixation3$choice == "correct") %>%
    ungroup() %>% complete(Condition, fix, choice) %>%
    ggplot() + geom_point(position = position_dodge(width = .8), aes(x = fix, y = n, color = Condition)) +
    facet_grid(Conf ~ .) +
    stat_summary(fun.y = "mean", geom = "crossbar", position = position_dodge(width = .8), 
                 mapping = aes(x = fix, y = n, color = Condition)) + ggtitle("fixation frequency in correct trials")

Stimulus-based fixation frequency (choice considered)

dat %>%
    group_by(Condition, fix, choice, id) %>%
    summarise(n = n()) -> fixation2
## `summarise()` has grouped output by 'Condition', 'fix', 'choice'. You can
## override using the `.groups` argument.
fixation2$fix <- as.factor(fixation2$fix)
fixation2$fix <- factor(fixation2$fix, levels = c("target", "distractor", "dud", "noFix"))

# correct choice
subset(fixation2, fixation2$choice == "correct") %>%
    ungroup() %>% complete(Condition, fix, choice) %>%
    ggplot() + geom_point(position = position_dodge(width = .8), aes(x = fix, y = n, color = Condition)) +
    stat_summary(fun.y = "mean", geom = "crossbar", position = position_dodge(width = .8), 
                 mapping = aes(x = fix, y = n, color = Condition)) + ggtitle("fixation frequency in correct trials")

# distractor choice
subset(fixation2, fixation2$choice == "distractor") %>%
    ungroup() %>% complete(Condition, fix, choice) %>%
    ggplot() + geom_point(position = position_dodge(width = .8), aes(x = fix, y = n, color = Condition)) +
    stat_summary(fun.y = "mean", geom = "crossbar", position = position_dodge(width = .8), 
                 mapping = aes(x = fix, y = n, color = Condition)) + ggtitle("fixation frequency in distractor-chosen trials")

# dud choice
subset(fixation2, fixation2$choice == "dud") %>%
    ungroup() %>% complete(Condition, fix, choice) %>%
    ggplot() + geom_point(position = position_dodge(width = .8), aes(x = fix, y = n, color = Condition)) +
    stat_summary(fun.y = "mean", geom = "crossbar", position = position_dodge(width = .8), 
                 mapping = aes(x = fix, y = n, color = Condition)) + ggtitle("fixation frequency in dud-chosen trials")
## Warning: Removed 6 rows containing non-finite values (`stat_summary()`).
## Warning: Removed 6 rows containing missing values (`geom_point()`).

Stimulus-based fixation proportion

fixation_prop1 <- c()

for (i in unique(fixation$id)) {
    for (j in unique(fixation$Condition)) {
        d <- subset(fixation,  fixation$id == i & fixation$Condition == j)
        s <- sum(d$n)
        fixation_prop1 <- rbind(fixation_prop1, mutate(d, prop = n/s))
    }
}

ggplot(fixation_prop1, aes(x = fix, y = prop, color = Condition)) + geom_violin() + 
    geom_point(position = position_dodge(width = .9)) +
    stat_summary(fun.y = "mean", geom = "crossbar", position = position_dodge(width = .9)) 

# noFix excluded)
fixation2 <- subset(fixation, fixation$fix != "noFix")
fixation_prop2 <- c()

for (i in unique(fixation2$id)) {
    for (j in unique(fixation2$Condition)) {
        d <- subset(fixation2,  fixation2$id == i & fixation2$Condition == j)
        s <- sum(d$n)
        fixation_prop2 <- rbind(fixation_prop2, mutate(d, prop = n/s))
    }
}

ggplot(fixation_prop2, aes(x = fix, y = prop, color = Condition)) + geom_violin() + 
    geom_point(position = position_dodge(width = .9)) +
    stat_summary(fun.y = "mean", geom = "crossbar", position = position_dodge(width = .9)) 

# dud excluded
fixation3 <- subset(fixation, fixation$fix == "target" | fixation$fix == "distractor")
fixation_prop3 <- c()

for (i in unique(fixation3$id)) {
    for (j in unique(fixation3$Condition)) {
        d <- subset(fixation3,  fixation3$id == i & fixation3$Condition == j)
        s <- sum(d$n)
        fixation_prop3 <- rbind(fixation_prop3, mutate(d, prop = n/s))
    }
}

ggplot(fixation_prop3, aes(x = fix, y = prop, color = Condition)) + geom_violin() + 
    geom_point(position = position_dodge(width = .9)) +
    stat_summary(fun.y = "mean", geom = "crossbar", position = position_dodge(width = .9)) 

fixation dynamics

df <- foreach(i = unique(dat$id), .packages = "tidyverse") %dopar% {
    df1 <- c()
    df2 <- subset(dat, dat$id == i)
    for (j in unique(df2$trial)) {
        df1 <- bind_rows(df1, df2 %>% filter(trial == j) %>% mutate(, nFix = row_number()))
    }
    print(df1)
}

dat2 <- c()
for (i in unique(dat$id)) {
    dat2 <- rbind(dat2, df[[i]])
}


dat2 %>%
    group_by(id, fix, nFix) %>%
    summarise(n = n()) -> fixs
## `summarise()` has grouped output by 'id', 'fix'. You can override using the
## `.groups` argument.
fix_prop <- c()
for (i in unique(fixs$id)) {
    for (j in 1:20) {
        fixs %>% filter(, id == i & nFix == j) -> d
        s <- sum(d$n)
        fix_prop <- rbind(fix_prop, mutate(d, prop = n/s))
    }
}

ggplot(fix_prop, aes(x = nFix, y = prop, color = fix)) + geom_point(position = position_dodge(width = .9)) +
    stat_summary(fun.y = "mean", geom = "line", position = position_dodge(width = .9))

fix_prop2 <- c()
for (i in unique(fixs$id)) {
    for (j in 1:20) {
        fixs %>% filter(, id == i & nFix == j & fix != "noFix") -> d
        s <- sum(d$n)
        fix_prop2 <- rbind(fix_prop2, mutate(d, prop = n/s))
    }
}

ggplot(fix_prop2, aes(x = nFix, y = prop, color = fix)) + geom_point(position = position_dodge(width = .9)) +
    stat_summary(fun.y = "mean", geom = "line", position = position_dodge(width = .9))